Groundwater level prediction based on a combined intelligence method for the Sifangbei landslide in the Three Gorges Reservoir Area

The monitoring and prediction of the groundwater level (GWL) significantly influence the landslide kinematics. Based on the long-term fluctuation characteristics of the GWL and the time lag of triggering factors, a dynamic prediction model of the GWL based on the Maximum information coefficient (MIC) algorithm and the long-term short-term memory (LSTM) model was proposed. The Sifangbei landslide in the Three Gorges Reservoir area (TGRA) in China, wherein eight GWL monitoring sensors were installed in different locations, was taken as a case study. The monitoring data represented that the fluctuation of the GWL has a specific time lag concerning the accumulated rainfall (AR) and the reservoir water level (RWL). In addition, there were spatial differences in the fluctuation of the GWL, which was controlled by the elevation and the micro landform. From January 19, 2015, to March 6, 2017, the measured data were used to set up the predicted models. The MIC algorithm was adopted to calculate the lag time of the GWL, the RWL, and the AR. The LSTM model is a time series prediction algorithm that can transmit historical information. The Gray wolf optimization (GWO) algorithm was used to seek the most suitable hyperparameter of the LSTM model under the specific prediction conditions. The single-factor GWO-LSTM model without considering triggering factors and the support vector machine regression (SVR) model were considered to compare the prediction results. The results indicate that the MIC-GWO-LSTM model reached the highest accuracy and improved the prediction accuracy by considering the factor selection process with the learner training process. The proposed MIC-GWO-LSTM model combines the advantages of each algorithm and effectively constructs the response relationship between the GWL fluctuation and triggering factors; it also provides a new exploration for the GWL prediction, monitoring, and early warning system in the TGRA.

The change of hydraulic condition is one of the crucial factors of landslides. As the most critical hydropower facility of the Yangtze River in China, the Three Gorges reservoir has dramatically changed the geological environment of this area 1 . The length of the TGRA is 637 km; therein more than 2500 landslides activated due to the fluctuation of the RWL 2 . This latter external factor, combined with the seasonal rainfall infiltration, significantly influences the seepage characteristics of landslides and the GWL 3 . The rise of the GWL increases the pore water pressure between the filled soil particles and the rock fissures and decreases the effective stress in the deformation area. The rapid decline of the GWL increases the dynamic water pressure, which accelerates the development of slope instability and deformation. The reduction in effective stresses due to the fluctuation of the GWL, affecting the soil strength and consequently the stability of the slope, causes first-time failures or the reactivation of landside movements 4 . Therefore, the monitoring and prediction of the GWL is an essential part of the landslide hazard analysis 5 . There is a correlation between the landslide displacement velocity and the GWL 6,7 . Van Asch 8 considered that the increase of the GWL could accelerate the deformation of landslides more than the decrease effect. On the contrary, Zhang 9 regarded the landslide deformation that occurs in the rapid GWL falling in the TGRA. The GWL changes the material strength at different locations of the landslide, which certainly accelerate Materials Geological and kinematic features. The Sifangbei landslide is located in the Wanzhou District, Chongqing (108° 29′ 18.65″ E, 30° 51′ 45.69″ N), which is on the north bank of the Yangtze River (Fig. 1). The Wanzhou District has abundant rainfall with a typical subtropical humid monsoon climate. According to the monitoring data of the Geo-Environmental Monitoring Office, the annual average temperature in Wanzhou District is 18.1 °C, the minimum temperature is − 3.7 °C, the maximum temperature is 42.1 °C, the relative humidity is 80%, the frost-free period is 334 days, and the annual sunshine hours are 1484.4 h. The average annual rainfall in Wanzhou is 1202. 8 mm (1960-2015). The rainfall mainly concentrates from May to September every year, accounting for about 70% of the annual rainfall (Fig. 2). The maximum average monthly rainfall is 203.1 mm in July. The historical maximum annual rainfall was 1635.2 mm (1981), the maximum monthly rainfall was 711.8 mm (July 1982), and the maximum daily rainfall was 243.31 mm (July 16,1982).
The Sifangbei landslide has an armchair shape in plain view with an elevation between 110 and 325 m. It has an estimated volume of 9.03 × 10 6 m 3 and covers an area of 3.612 × 10 5 m 2 , with a length of 840 m and a width of 430 m. The sliding body depth (Fig. 2) varies from 10.3 to 28.9 m, and the average bedrock depth is 22 m. The surface water network of Sifangbei landslide is relatively developed, mainly including a water pond, farmland, and gully. The pond on the landslide has water all year round, with a total area of about 1.2 × 10 4 m 2 . According to the occurrence characteristics, the groundwater of the Sifangbei landslide can be divided into bedrock fissure water and loose pore water. The pore water is mainly distributed in the Quaternary deposits and rubble soil layer, mainly supplied by atmospheric precipitation and surface water, and discharged into the river along the bedrock surface.
The engineering geological profile of the Sifangbei landslide is stepped form with a main sliding direction in 160° (Fig. 3). There are two-level platforms developed in the landslide area, the central platform is 220 m long, 450 m wide, and 200 m high. The sliding surface in correspondence of the landslide toe is below the water level of the Yangtze River. The surface layer of the landslide consists of thick quaternary silty clay mixed with broken rocks. The underlying layer is typical near-horizontal Jurassic strata in Wanzhou District. The occurrence of the strata is 160°∠5°.
The initial monitoring network of the Sifangbei landslide was established in March 2007, which includes 3 GPS stations (Fig. 1) www.nature.com/scientificreports/ the Sifangbei landslide was accelerated due to the first water storage of Three Gorges Reservoir to 156 m, which caused many through cracks and secondary collapses in the middle-front area of the landslide. Based on the surface displacement data collected from March 2007 to April 2016 (Fig. 4), the cumulative displacement of GPS03 is the largest, followed by GPS02 and GPS01. The rapid increase of displacement mainly occurred during the rapid   www.nature.com/scientificreports/ decline period of RWL and the month of heavy rainfall. Accordingly, it can be concluded that the deformation of the Sifangbei landslide is mainly affected by RWL and rainfall.
In-situ monitoring and data acquisition. Members of the research group started a field investigation in Sifangbei in July 2014 and established a hydrological monitoring system in 2015. Data for the present study are collected from the Sifangbei landslide hydrological monitoring system. The monitoring system (see Fig. 1) consists of six Semi-automatic GWL holes (STK1-STK6) and two automatic GWL holes (STK7, STK8). The semiautomatic GWL holes are measured by a MicroDiver probe, which can store up to 48,000 groups of data, with an accuracy of ± 0.05%, and a regular working temperature of 0-40 °C. The automatic GWL holes are a voltage-type water level gauge, with an accuracy of ± 0.25%. Both holes measure GWL once a day. After collecting and saving the device record data, the GWL can be obtained through the calculation formula: where H is the groundwater level (m), h is the device record data (cm), h s is the height of the water column at standard atmospheric pressure (m) and it is equal to 10.336 m, H s is the elevation of the monitoring holes (m), d is the distance between the MicroDiver probe and the monitoring hole.
(1)  www.nature.com/scientificreports/ The GWL of the Sifangbei landslide has been monitored since 2015, including the STK1-STK6 from January 19, 2015, and the STK7, and STK8 from October 30, 2018. Considering the missing data and maintenance of monitoring points, periods with more complete data were selected for analysis. The data of the semi-automatic GWL holes were selected from January 19, 2015, to March 6, 2017, and data of the automatic GWL holes were selected from October 30, 2018, to May 28, 2020 (Fig. 5).
The GWL monitoring holes can be divided into four groups ( Table 1). The distance of each group from the Yangtze River is shown in Fig. 1. According to different groups, the GWL of Sifangbei can be analyzed: (1) The STK7 and STK8 are located in the fluctuation range of RWL (175-145 m). The fluctuation of monitoring data of these two points is positively correlated with the RWL, and almost at the same value. The monitoring point is the underwater part of each year, so it is less affected by rainfall. (2) The elevation of the STK1 and STK4 is 175 m and they are close to the Yangtze River. The fluctuation of the GWL is affected by the RWL. The GWL of STK1 ranges from 169.69 to 173.59 m as a gentle slope behind the monitoring hole (Fig. 6a) and also shows a short-term fluctuation during the rainy season. The GWL of STK4 ranges from 168.29 to 177.88 m, mainly due to the loose surface soil and high permeability coef- www.nature.com/scientificreports/ ficient caused by the planting of fruit trees (Fig. 6b). The increase of hydrodynamic pressure and seepage pressure accelerates the deformation of the landslide front area. (3) The STK2 and STK5 are located on the first-level platform of the Sifangbei landslide. The STK2 is located near farmlands and ponds and is mainly affected by human agricultural activities. The GWL of STK2 is always maintained at about 180 m in the whole process of RWL regulation. The GWL of STK5, mainly affected by rainfall, is far away from the Yangtze River and mainly fluctuates in the rainy season. (4) The STK3 and STK6 are the farthest monitoring points from the Yangtze River and are mainly controlled by rainfall. The platform on the slope collected water in the rainy season, resulting in the rapid rise of the GWL. After the rainy season, the groundwater was discharged along the slope. The primary performance of these two points represents a periodic fluctuation within the year, which is a rapid rise followed by a rapid fall.
The GPS03, STK5, and STK6 hydrological holes (relatively close) were selected to analyze the relationship between the horizontal surface displacement and the GWL. As shown in Fig. 7, there is a relatively good fitting between the surface displacement and the GWL. The fluctuation of the GWL causes the fluctuation of cumulative displacements. The adjacent elevation hydrological hole (STK5) can better reflect the fluctuation in surface displacement (GPS03). The rapid fluctuation of the GWL causes an increase in displacement. Especially from February 2015 to September 2015, the cumulative displacement increased by 77.1 mm. Therefore, the monitoring and analysis of groundwater levels can reflect the characteristics of landslide deformation. Table 1 represents the qualitative correlation of different monitoring points with rainfall and RWL. At the front edge of the landslide, the fluctuation of the GWL is mainly affected by RWL, and at a distance from the Yangtze River, it is mainly affected by rainfall. The STK1 monitoring hole, which the RWL and rainfall influence, was used to establish the prediction model for the present study.

Methodology
Time series analysis. In the TGRA, the GWL in reservoir landslide is related to the RWL and the rainfall.
According to the data of the RWL, the AR, and the GWL, the time-series relationship is established as follows:  www.nature.com/scientificreports/ where g t is the GWL at time t, f is the prediction model, x n is RWL factors, y i is rainfall factors.
Long short-term memory neural networks. The input-implicit-output layer is fully connected to traditional neural networks. The points between the sequences are not connected. Therefore, the traditional neural network is not suitable for time series prediction 48 . The recurrent neural network (RNN) records the previous information and applies it to the current output calculation, which allows for feedback in the networks 49 . The LSTM neural networks are a special type of RNN, which can effectively solve the problem of "gradient disappearance" or "gradient explosion" 50 . The LSTM ensures data discovery and long-term memory through the forgetting gate, input gate, and output gate. The three gate functions provide a good nonlinear control mechanism for the input and deletion of control information (Fig. 8). The hidden vector (h (t) ) in the LSTM model can be obtained as follows: where f (t) , i (t) , o (t) and c '(t) are the values of the forget gate, the input gate, forget gate, the output gate, and the memory cell in the memory block; and b o are their corresponding bias values; σ is the sigmoid function, tanh is the hyperbolic tangent function; C (t) is the updated value of the cell state; y (t) is the output value at time t. After forward propagation, BPTT (backpropagation through time) algorithm is used to transfer the accumulated error back from the last time, calculate the gradient of error corresponding parameters. Finally, the weights and thresholds are updated by the stochastic gradient descent algorithm.
(3)  www.nature.com/scientificreports/ Grey wolf optimizer. The GWO is a novel meta-heuristic algorithm with strong search ability that imitates the leadership and hunting mechanism of gray wolves in nature 51 . It uses three wolf heads (α wolf, β wolf, and δ Wolf) to determine the fitness value, while other wolves calculate the distance between them and their prey according to the location of their prey 51 (Fig. 9).
(1) According to the level of the wolf group, the optimal solution is regarded as α wolf, the second and third optimal solutions are β wolf and δ wolf respectively, and the other candidate solutions are ω wolf.
(2) The distance D between the prey and the grey wolves is determined before preying: where C is the swing factor, � C=2� r ; r is the random vector,� r = random[0, 1][0, 1] ; X P (t) is the position vector of the prey of the t th generation of grey wolves; X(t) is the position vector of the t th generation of grey wolves.
(3) The distance between the gray wolf and its prey is shortened by iterative updating: where A is the convergence factor; a is the weight factor, and the initial value is 2, which decreases to 0 as the number of iterations increases. A grey wolf in the position of (X, Y) can update its position according to the position of the prey (X*, Y*). Different places around the best agent can be reached concerning the current position by adjusting the value of A and C vectors.
(4) In the abstract search space, the ω wolf moves closer to the prey according to the position of the three wolves (Fig. 10). The final position is determined by the random position in the circle defined by α, β, and δ in the search space. The distance between the other grey wolves and these three wolves in the t the generation can be obtained according to the following formula:  www.nature.com/scientificreports/ The number of hidden layer neurons, the number of iterations, and the learning rate have a significant impact on the fitting and prediction ability of the LSTM model. Therefore, the GWO algorithm was introduced to optimize the hyperparameter of the LSTM model in landslide displacement prediction. The flowchart is depicted in Fig. 11. The optimal hyperparameters of the GWO-LSTM model were obtained by using the training set and the testing set. Considering the continuity and historical memory of time series, the training set and the testing set were combined as the new training set, and they get the prediction results compared with the verification set, which verifies the prediction accuracy and generalization ability of the model.

Maximal information coefficient. The MIC is a new variable correlation analysis method proposed by D.
n. Reshef et al. 52 in 2011 based on mutual information (MI). This method can not only find the linear functional relationship between variables but also find the nonlinear functional relationship (exponential, periodic, etc.) 44 . The main steps of the algorithm are as follows: (1) Supposed there is a connection between two variables V 1 = {v 1 (i)}, i = 1,2,…,n and V 2 = = {v 2 (i)},i = 1,2,…,n.
D is the set of ordered pairs {v 1 (i),v 2 (i)},i = 1,2,…,n. Using grid G 1 (x 1 × y 1 ) to divide V 1 sample points into x 1 and V 2 sample points into y 1 , some cells are allowed to be empty sets. (2) The characteristic matrix M(D) x 1 ,y 1 can be obtained from the maximum mutual information value max I(D| G 1 ) as follows: where D| G 1 is the probability mass distribution function of all cells in grid G 1 ; n ij is the sample point in column i of row j in grid G 1 , N is the total number of samples. Because different grid G will lead to different D| G , the global optimal grid G 0 (x 0 × y 0 ) is determined by an exhaustive search of the characteristic matrix. The maximum information coefficients of variables V 1 and V 2 are as follows: where B (N) is the maximum grid area for searching. In general, The MIC algorithm is essentially a normalized maximum mutual information, and its value range is [0, 1]. For two independent variables, the MIC tends to be 0; for two variables with noise-free function relationships, the MIC tends to 1. To verify the accuracy of prediction, the RMSE, and the R 2 was used to evaluate the prediction accuracy 53 .  www.nature.com/scientificreports/ ber-December). Accordingly, the hydrogeological conditions of the TGRA present different characteristics over the year 30 . The long-term fluctuation of the GWL is related to the RWL, while the short-term fluctuation is affected by AR. The GWL represents different degrees of time lag with the change of external conditions and soil permeability. It can be concluded that the response time lag of the STK1 monitoring GWL to the RWL is short from the geographical location and monitoring data.

Results
The most crucial step of the GWL prediction is determining the relevant triggering factors. The MIC algorithm was adopted to determine the correlation between the GWL and triggering factors. Figure 13a shows that the correlation is the largest when the RWL leads GWL 28-35 days. The curve chart of the RWL and the GWL (Fig. 14) shows that the overall trend is consistent with the leading RWL; the peak value of the GWL is the same as that of the current RWL, the valley value periodic change is consistent with that of the leading RWL.
The short-term fluctuation of the GWL is affected by the AR. The AR and the GWL data from April 1, 2016, to September 1, 2016, were selected for correlation analysis (Fig. 13b). In particular, the most significant correlation  www.nature.com/scientificreports/ was recorded for 9 and 13 AR days with values of 0.556 and 0.563, respectively. As shown in Fig. 15, the GWL raised with the increase of AR. There is a strong correlation between the GWL and the AR. In summary, the RWL leading 28 days and 35 days, and the AR leading 9 days and 13 days were chosen as the optimal triggering factors of the GWO-LSTM model.
Prediction results and analysis. Seven hundred and seventy-eight measured data from January 19, 2015, to March 6, 2017, were adopted to establish the predicted model. Six hundred and twenty-two data from January 19, 2015 to October 1, 2016 were selected as the training samples. The following 78 groups of data (October 2, 2016-December 18, 2016) were used as the test samples to evaluate the model's accuracy and determine the optimal prediction model. The data from December 19, 2016 to March 6, 2017 were used as a validation set to verify the generalization ability of the optimal prediction model. The optimal hyperparameters of the LSTM model were obtained by using the training set and the testing set. Then, the training set and the testing set were combined into the training set of the validation set. The maximum and minimum values of input and output for normalization and inverse normalization were obtained respectively from the training data set to prevent information leakage. Each triggering factor was normalized to [− 1, 1] using linear normalization: where x * is the normalized value, x is the original value, x max is the maximum value of the samples, and x min is the minimum value of the samples. Based on an original program written in Matlab R2020b for this specific analysis, the LSTM model and the SVR model were used to learn and train the data in the training samples. When the network error converged to the expected value, the network was used to predict the data in the test samples to test the generalization ability. The single-factor GWO-LSTM model, the MIC-GWO-LSTM model, and the MIC-GWO-SVR model were proposed in this paper to compare and analyze the prediction performance of different models.
As  36 . In the GWO algorithm parameter setting, the penalty factor C = [0,100], the kernel function parameter is γ = [0, 100], the principal component is set to 95%, and the cross-validation value is v = 5. The number of grey wolf groups and generations were set as 30, 50 respectively. The RMSE between the measured and predicted values was set as the fitness function.
Prediction results of the test set, using the GWO algorithm to search for optimal hyperparameters, are shown in Table 2. The three optimized models show good prediction ability on the test set (Fig. 16a). The RMSE and R 2 of the single-factor GWO-LSTM model were 0.0624 and 0.9974. The RMSE and R 2 of the multi-factor GWO-SVR model were 0.0615 and 0.9975, whereas the values of the MIC-GWO-LSTM model were 0.0544 and 0.9977, respectively. The MIC-GWO-LSTM model has a better training ability than the MIC-GWO-SVR and the single-factor GWO-LSTM model. The most significant errors for the single-factor GWO-LSTM (Relative error 0.125%), MIC-GWO-SVR (Relative error 0.187%), and MIC-GWO-LSTM (Relative error 0.096%) models for the direct prediction of the GWL occurred on November 8 (Fig. 16b,c). This time point is the inflection point when the RWL rose from 145 to 175 m. It shows that these prediction models have a certain prediction delay for the inflection point of data.
The optimized hyperparameters (Table 2) were used to predict the validation set. Figure 17a shows that the predicted results of the three models were all consistent with the trend of measured data. Due to the lack of triggering factors, the single-factor GWO-LSTM model, which represents more fluctuation in the validation set, has poor prediction and generalization ability (RMSE = 0.0660, R 2 = 0.9641). The RMSE and R 2 of the MIC-GWO-SVR model were 0.0532 and 0.9769, whereas the values of the MIC-GWO-LSTM model were 0.0457 and 0.9830, respectively. The single-factor GWO-LSTM model has higher fluctuations and errors in the whole prediction sequence (Fig. 17b), and has a lower generalization ability than the MIC-GWO-LSTM model. The MIC-GWO-SVR model and the MIC-GWO-LSTM model have different prediction abilities and represent different prediction errors on the validation set (Fig. 17c). They represent the complicated nonlinear relationship between the GWL and its triggering factors and have good generalization ability on the validation set. In general, the MIC-GWO-LSTM model represents a higher prediction accuracy than the MIC-GWO-SVR model.  www.nature.com/scientificreports/

Discussion
The accurate prediction of the GWL needs to find the main control factors 55 . The AR and the RWL are the two main factors that can change the hydrological conditions of the GWL in the TGRA. Four sets of prediction experiments were proposed with different triggering factors to compare the MIC-based feature engineering's prediction stability. As shown in Table 3, models of four different triggering factors were established and compared with the MIC-GWO-LSTM model. They all have good prediction ability in the training set, especially in the stage of groundwater level fluctuation from October 27 to November 14 (Fig. 18), mainly as the GWO algorithm searched the suitable hyperparameter of the LSTM model under specific prediction conditions. The prediction results of the four models are quite different in the validation set (Fig. 19). The prediction accuracy and generalization ability of the MIC-GWO-LSTM model decreased with input factors. These triggering factors contain a lot of duplicates and redundant information because they are mainly calculated by the measured data of the current  www.nature.com/scientificreports/ day. Two factors (model-4) with high correlation were screened out by the Pearson correlation coefficient (PCC) algorithm. The model-4 also predict the GWL with high prediction and generalization ability. The redundant triggering factors will reduce the prediction accuracy. The length of input days is not the only criterion for input factors, and the interaction of different factors should also be considered. Krakac 4,26 selected many triggering factors to predict groundwater depth, including 0-100-day antecedent precipitation and 0-100-day evapotranspiration. This method only focuses on the data itself and ignores the natural phenomena. More input data will increase the operation time and not guarantee prediction accuracy. Cao 12 applied the grey correlation model to calculate the correlation degree between groundwater and triggering factors but ignored the collinearity of redundant features, which may reduce the generalization ability of the prediction model. In this paper, the time lag of the GWL relative to the RWL and the RA was determined by the MIC-based feature engineering, which integrated the factor selection process with the learner training process, and effectively improved the prediction accuracy and the generalization ability of the prediction model. One should consider developing accurate triggering factors for reducing the model uncertainly in the future. The LSTM model saves and uses the historical information and gives full play to extract correlation information 45 . It also represents a good prediction result in the situation of enough samples. The main control factor of the LSTM model is the hyperparameters. Compared with the grid search, the GA, and the PSO algorithms, the GWO algorithm avoids falling into local optima in high-dimensional and has a lower convergence Table 3. The prediction accuracy comparisons between different triggering factors of the GWO-LSTM model.

Model
The triggering factors set www.nature.com/scientificreports/ speed in the iterative process. In this paper, The GWO-LSTM model has better generalization performance and improves the prediction accuracy of the GWL by searching for the hyperparameters of the LSTM model. The fluctuation of the GWL is a complex dynamic system, which is related to hydrological and geological conditions. In the long-term observation period, the effect of periodicity and randomness is reflected in the GWL measured data. The GWL is also closely related to the slope structure and macroscopic deformation characteristics (earth cracks, etc.). Earth cracks increase the permeability of the soil and make groundwater flow more frequently 12 . How to predict the GWL by combining the characteristics of meteorology, hydrology, landslide mass, and landslide macro deformation needs further research.
Compared with the RA and the RWL, the GWL can directly reflect the hydrological characteristics of landslides. The landslide velocity changed with the fluctuation of the GWL because it causes the change of pore water pressure in the soil and directly affects the deformation characteristics 10 . Therefore, a more practical GWL prediction is helpful to predict landslide displacements, which could be valuably exploited to establish a landslide early warning system.

Conclusion
The monitoring and prediction of the GWL are essential for establishing a landslide early warning system. A MIC-GWO-LSTM model for predicting the GWL of the landslide was tested. The conclusions are summarized below: 1. There is a specific correlation in time series between the GWL with the RA and the RWL. The influence of RWL on the Sifangbei landslide is mainly concentrated on the front area. The GWL is mainly affected by rainfall as the distance from the Yangtze River increases. The farmland and pond also affect the change of the GWL. 2. The MIC algorithm effectively evaluates the relationship between factors and the GWL and analyses the main triggering factors of the RWL and the AR, which improved the generalization ability of prediction modeling. 3. The results show that the prediction accuracy of the GWL prediction model based on the MIC-GWO-LSTM model is better than that of other prediction models. www.nature.com/scientificreports/ In general, the MIC-GWO-LSTM model proposed in this paper can effectively construct the response relationship between the GWL and triggering factors. Meanwhile, many machine learning algorithms have been developed for landslide susceptibility mapping and risk analysis 27,56 . Thus, the novel method proposed in this paper is recommended for conducting landslide susceptibility mapping, landslide risk analysis, and other fields.

Data availability
Restrictions apply to the availability of these data. Data were obtained from the Geo-environmental monitoring office of Wanzhou District and are available from the authors with the permission of the Geo-environmental monitoring office of Wanzhou District.